The purpose of the first set of exercises is, on one hand, to familiarize with the programming environment of Python and Matlab, and on the other hand, to introduce the methods of representation and processing of telecommunication signals in these specific programming languages.
Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).
fromscipyimportsignal# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal# https://docs.scipy.org/doc/scipy/reference/signal.htmlfromscipy.fftimportfft,fftfreqfromscipy.fftpackimportfftshift,ifftshiftimportmatplotlibimportmatplotlib.pyplotaspltimportnumpyasnpfromnumpyimportrandomimportwarningswarnings.filterwarnings('ignore')fromipywidgetsimportIntRangeSlider,widgets,Layout,HBox,VBoxfromIPython.displayimportdisplay,clear_outputprint("Libraries added successfully!")
a[0:2,:]# Ομοίως: οι γραμμές 1 & 2 δίνονται ως 0:2 και όχι ως 0:1
a_slice=a(1:2,:);disp(a_slice);
Create a vector with elements from 0 to 0.5 with a step of 0.1#
t=np.arange(0,0.5,0.1)print('t=',t)
t=0:0.1:0.5;% The endpoint is inclusive in MATLAB, so you might want to adjust it if you want the exact same output as NumPyt(end)=[];% Remove the last element to match Python's exclusive endpointdisp(['t = ',mat2str(t)]);
Primary signals are mainly analog (continuous time). To represent and process them on our computer (or another digital machine), we must first digitize them. Consider a continuous-time signal \(x(t)\) with Fourier transform (Continuous Time Fourier Transform – CTFT):
# Parametersfs=1000# Sampling frequency in HzT=1/fs# Sampling periodL=1000# Number of samplest=np.arange(0,L)*T# Create time vector# Output widget for the plotsgraph_output0=widgets.Output()defupdate_plot(selected_frequencies_first):withgraph_output0:clear_output(wait=True)# Clear the previous plot# Generate the signal with the selected frequenciessignal=np.zeros(len(t))forfinrange(selected_frequencies_first[0],selected_frequencies_first[1]+1):signal+=np.sin(2*np.pi*f*t)# Calculate the FFT of the original signaloriginal_signal_fft=fft(signal)# Calculate frequencies for the FFT of the original signaloriginal_freqs=fftfreq(L,T)[:L//2]# Calculate magnitude of Fourier coefficients (amplitude) for the original signaloriginal_magnitude=np.abs(original_signal_fft)[:L//2]# Plottingfig,axs=plt.subplots(1,2,figsize=(15,5))# Original signal plotaxs[0].plot(t,signal,color='#00CC96')axs[0].set_title('Original Signal')axs[0].set_xlabel('Time (s)')axs[0].set_ylabel('Amplitude')axs[0].grid(True)# Fourier Transform plotaxs[1].plot(original_freqs,original_magnitude)axs[1].set_title('Fourier Transformation of Original Signal')axs[1].set_xlabel('Frequency (Hz)')axs[1].set_ylabel('Magnitude')axs[1].grid(True)plt.tight_layout()plt.show()# Create the slider widgetfrequency_slider=widgets.IntRangeSlider(value=[5,10],min=1,max=50,step=1,description='Frequency Range (Hz)',style={'description_width':'initial'},layout=Layout(width='90%'),continuous_update=False)defresponse(change):fig=update_plot(frequency_slider.value)# Observe changes in the slider and update the plot accordinglyfrequency_slider.observe(response,names='value')# Style the HTML elementhtml_label=widgets.HTML(value=""" <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Continuous Time Fourier Transform – CTFT</h2> """)# Display the slider and the outputsvbox_layout=Layout(display='flex',flex_flow='column',align_items='center',width='100%')ui=widgets.VBox([html_label,frequency_slider,graph_output0],layout=vbox_layout)out=widgets.interactive_output(update_plot,{'selected_frequencies_first':frequency_slider})clear_output(wait=True)# Clear the previous plotdisplay(ui,out)
By sampling \(x(t)\) at a rate \(f_s=1/T_s\), a discrete-time signal \(x(nT_s)\) is produced. Mathematically, it is represented as a series of delta functions
# Παράμετροι δειγματοληψίαςT=1/fs# Sampling periodL=1000# Number of samples# Create time vectort=np.arange(0,L)*Tdefplot_sampled_signal(f_range,n_samples):# Generate the signal with the selected frequenciessignal=np.zeros(len(t))forfinrange(f_range[0],f_range[1]+1):signal+=np.sin(2*np.pi*f*t)# Correct downsampling approachsample_rate=n_samples# Directly using selected_samples as an integerdownsampling_factor=int(fs/sample_rate)# Downsampling factor# Select every nth sample from the original time vector to match the downsampled signalsample_points=t[::downsampling_factor]# Downsample the signal by selecting every nth sampledownsampled_signal=signal[::downsampling_factor]# Recalculate L for the downsampled signal if neededL_downsampled=len(downsampled_signal)# Calculate the DFT of the downsampled signalsampled_signal_fft=fft(downsampled_signal)# Calculate frequencies for the FFT of the downsampled signal# Note: The new sampling period is the inverse of the new sampling ratesampled_freqs=fftfreq(L_downsampled,1/sample_rate)[:L_downsampled//2]# Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signalsampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]n=len(downsampled_signal)T=1/sample_rate# Recalculate the sampling period with the selected sample rate# After calculating the magnitude of Fourier coefficientssampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]# Create a new frequency vector that includes the replicated frequencies# First, calculate the original frequency bins for the positive frequenciesoriginal_sampled_freqs=fftfreq(L_downsampled,T)[:L_downsampled//2]shifted_freqs_negative=original_sampled_freqs-1/T# Shift for -1/Tsshifted_freqs_positive=original_sampled_freqs+1/T# Shift for +1/Ts# Concatenate the original and shifted (replicated) frequencies and magnitudes# This includes the original frequencies, and the replications at -1/Ts and +1/Tsreplicated_freqs=np.concatenate([shifted_freqs_negative,original_sampled_freqs,shifted_freqs_positive])replicated_magnitude=np.concatenate([sampled_magnitude,sampled_magnitude,sampled_magnitude])# Sort the replicated frequencies and magnitudes in ascending order for proper plottingsorted_indices=np.argsort(replicated_freqs)sorted_replicated_freqs=replicated_freqs[sorted_indices]sorted_replicated_magnitude=replicated_magnitude[sorted_indices]# Plottingfig,axs=plt.subplots(1,2,figsize=(15,5))axs[0].scatter(sample_points,downsampled_signal,color='#00CC96',s=12)axs[0].set_title('Sampled Signal')axs[0].set_xlabel('Time (s)')axs[0].set_ylabel('Amplitude')axs[0].grid(True)axs[1].plot(sorted_replicated_freqs,sorted_replicated_magnitude,color='#1F77B4')axs[1].set_title('Fourier Transformation of Sampled Signal')axs[1].set_xlabel('Frequency (Hz)')axs[1].set_ylabel('Magnitude')axs[1].grid(True)plt.tight_layout()plt.show()# Widgetsfrequency_slider=widgets.IntRangeSlider(value=[5,10],min=1,max=50,step=1,description='Frequency Range (Hz):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)samples_slider=widgets.IntSlider(value=500,min=100,max=1000,step=25,description='Samples:',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Signal Reconstruction in Digital Signal Processing</h2> """)vbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,frequency_slider,samples_slider],layout=vbox_layout)out=widgets.interactive_output(plot_sampled_signal,{'f_range':frequency_slider,'n_samples':samples_slider})clear_output(wait=True)# Clear the previous plotdisplay(ui,out)
For band-limited signals \(x(t)\) with bandwidth W, assuming that the sampling rate \(fs ≥ 2W\), it holds that \(X(f) = T_s X_\delta(f)\), \(0 ≤ f ≤ W\), meaning, the signal \(X(f)\) results from passing the sampled \(x_\delta(t)\) through an ideal low-pass filter with gain \(T_s\). From the previous schema, it becomes evident that if the sampling is done with a frequency less than twice the highest frequency \(W\) of the signal (undersampling), then “images” of the spectrum from higher frequencies appear in the signal’s frequency domain, preventing the accurate restoration of the original continuous-time signal. This phenomenon is called aliasing, and the error during the restoration of the original signal is referred to as aliasing error.
Sampling in the time domain is the basis for defining the Discrete Time Fourier Transform (DTFT). For a series of discrete numbers \(x[n]\), the discrete-time Fourier transform is defined as:
The DTFT is a periodic function with a period of \(1\), therefore, it is sufficient to compute it in the frequency interval \([0,1]\) or equivalently \([-½,½]\). Note that the DTFT, although it arises from a series of discrete numbers \(x[n]\), is a continuous function of the variable \(\phi\) as illustrated in the next schema.
# Παράμετροι δειγματοληψίαςT=1/fs# Sampling periodL=1000# Number of samples# Create time vectort=np.arange(0,L)*Tdefplot_sampled_signal(selected_frequencies,selected_samples):# Generate the signal with the selected frequenciessignal=np.zeros(len(t))forfinrange(selected_frequencies[0],selected_frequencies[1]+1):signal+=np.sin(2*np.pi*f*t)# Correct downsampling approachsample_rate=selected_samples# Directly using selected_samples as an integerdownsampling_factor=int(fs/sample_rate)# Downsampling factor# Downsample the signal by selecting every nth sampledownsampled_signal=signal[::downsampling_factor]# Recalculate L for the downsampled signal if neededL_downsampled=len(downsampled_signal)# Calculate the DFT of the downsampled signalsampled_signal_fft=fft(downsampled_signal)# Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signalsampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]n=len(downsampled_signal)T=1/sample_rate# Recalculate the sampling period with the selected sample rate# After calculating the magnitude of Fourier coefficientssampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]# Create a new frequency vector that includes the replicated frequencies# First, calculate the original frequency bins for the positive frequenciesoriginal_sampled_freqs=fftfreq(L_downsampled,T)[:L_downsampled//2]shifted_freqs_negative=original_sampled_freqs-1/T# Shift for -1/Tsshifted_freqs_positive=original_sampled_freqs+1/T# Shift for +1/Ts# Concatenate the original and shifted (replicated) frequencies and magnitudes# This includes the original frequencies, and the replications at -1/Ts and +1/Tsreplicated_freqs=np.concatenate([shifted_freqs_negative,original_sampled_freqs,shifted_freqs_positive])replicated_magnitude=np.concatenate([sampled_magnitude,sampled_magnitude,sampled_magnitude])# Sort the replicated frequencies and magnitudes in ascending order for proper plottingsorted_indices=np.argsort(replicated_freqs)sorted_replicated_freqs=replicated_freqs[sorted_indices]sorted_replicated_magnitude=replicated_magnitude[sorted_indices]# Generate an array of sample indices for the downsampled signalsample_indices=np.arange(len(downsampled_signal))# Normalize the frequency valuesnormalized_freqs=original_sampled_freqs/(1/T)# Since you previously concatenated and sorted for replication, ensure to apply normalization there as wellnormalized_replicated_freqs=np.concatenate([(shifted_freqs_negative/(1/T)),# Normalize -1/Ts shifted frequenciesnormalized_freqs,# Already normalized original frequencies(shifted_freqs_positive/(1/T))# Normalize +1/Ts shifted frequencies])# Sort the normalized and replicated frequencies for proper plottingsorted_indices=np.argsort(normalized_replicated_freqs)sorted_normalized_replicated_freqs=normalized_replicated_freqs[sorted_indices]sorted_replicated_magnitude=replicated_magnitude[sorted_indices]# Create the plotsfig,axs=plt.subplots(1,2,figsize=(15,5))# Plot the downsampled signalaxs[0].scatter(sample_indices,downsampled_signal,color='#00CC96',s=12)axs[0].set_title('Sampled Signal')axs[0].set_xlabel('Sample Number')axs[0].set_ylabel('Amplitude')axs[0].grid(True)# Plot the Fourier Transform of the downsampled signalaxs[1].plot(sorted_normalized_replicated_freqs,sorted_replicated_magnitude,color='#1F77B4')axs[1].set_title('Normalized Fourier Transformation of Sampled Signal')axs[1].set_xlabel('Normalized Frequency (f/fs)')axs[1].set_ylabel('Magnitude')axs[1].grid(True)plt.tight_layout()plt.show()# IPyWidgets slidersfrequency_slider=widgets.IntRangeSlider(value=[5,10],min=1,max=50,step=1,description='Frequency range (Hz):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)samples_slider=widgets.IntSlider(value=500,min=100,max=1000,step=25,description='Samples:',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Discrete Time Fourier Transform – DTFT</h2> """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,frequency_slider,samples_slider],layout=vbox_layout)# Interaction between widgets and functionout=widgets.interactive_output(plot_sampled_signal,{'selected_frequencies':frequency_slider,'selected_samples':samples_slider})clear_output(wait=True)# Clear the previous plotdisplay(ui,out)
With the series of discrete numbers arising as a result of sampling, \(x[n]=x(nT_s)\), the DTFT and the Fourier transform \(X_\delta(f)\) of the sampled signal are connected through the correspondence \(\phi ↔ f/f_s\). The common practice is to represent the ratio \(f/f_s\) as normalized frequency \(\phi\) (\(f_D\), in your notes) and the real frequencies to arise as multiples of it (usually fractional). To link the DTFT with the Fourier transform \(X(f)\) of the signal, an additional conversion to the sampling period by multiplying by \(T_s\) (or dividing by \(f_s\)) is required.
Analogous to sampling signals in time, we can sample in the frequency domain by taking discrete values \(X(kf_o)\) of the Fourier transform corresponding to a frequency analysis \(f_o=1/T_o\). This is equivalent to the periodic repetition of the continuous-time signal \(x(t)\) every \(Τ_ο\), since the periodic signal
Therefore, \(X[k] = X(kf_o)/Τ_o\) are the coefficients of the Fourier series expansion of the periodic signal \(x_p(t)\). Obviously, for signals \(x(t)\) of finite duration, where \(x(t)=0\) for \(|t| ≥ T\), assuming that the period \(T_o ≥ 2T\), it holds that \(x(t) = x_p(t)\) for \(|t| ≤ T\).
In practice, signals have a very long duration to be analyzed in their entirety. Thus, we apply a rectangular time window to retain only their most significant part for the observation period and \(x(t)= 0\) elsewhere. In computing the DTFT \(X_d(\phi)\) of such a truncated signal, instead of an infinite sum, we are limited to a finite-length \(L\) series of numbers \(x[n]\), therefore
\[
X_d(\phi)=\sum_{n=0}^{L-1} x[n] \exp (-j 2 \pi n \phi)
\]
The sampling of \(X_d(\phi)\) in the frequency domain at \(N\) equidistant normalized frequencies \(0\), \(1/N\), \(2/N\), \(…\), \((N-1)/N\), gives
\[
X[k]=X_{d}\left(\frac{k}{N}\right)=\sum_{n=0}^{N-1} x[n] \exp \left(-j 2 \pi n \frac{k}{N}\right), \quad 0 \leq k \leq N-1
\]
where, if \(N≥L\), we set \(x[n]=0\) for \(n≥L\). The last relationship is recognized as the Discrete Fourier Transform (DFT), which for a finite series \(xn\), \(n=0\), \(1\), \(…\), \(N-1\), is defined as:
\[
X_{k} \triangleq \sum_{n=0}^{N-1} x_{n} \exp \left(-j 2 \pi n \frac{k}{N}\right), \quad 0 \leq k \leq N-1
\]
and its inverse is
\[
x_{n}=\frac{1}{N} \sum_{k=0}^{N-1} X_{k} \exp \left(j 2 \pi n \frac{k}{N}\right), \quad 0 \leq n \leq N-1
\]
The \(X_d(\phi)\) as DTFT is a periodic function and if the original series xn was periodic (and we didn’t apply the window), then \(X_d(\phi)\) would be zero everywhere except at the sampling points \(k/N\). That is, if we consider a finite-length series of numbers that repeats periodically, the discrete-time Fourier transform of it (DTFT) is also periodic and discrete. Moreover, the DFT and its inverse IDFT, if we did not limit the indices \(n\) and \(k\) between \(0\) and \(N-1\), would be periodic functions. Therefore, the finite series xn can be considered as a periodic discrete-time signal observed only during a period and the DFT, the series \(X_k\), as the samples with resolution \(1/N\) of the DTFT \(X_d(\phi)\) in the domain of normalized frequencies \([0,1]\), as shown in the next schema.
# Παράμετροι δειγματοληψίαςT=1/fs# Sampling periodL=1000# Number of samples# Create time vectort=np.arange(0,L)*Tdefplot_signals(f_range,samples,selected_N):# Generate the signal with the selected frequenciessignal=np.zeros(len(t))forfinrange(f_range[0],f_range[1]+1):signal+=np.sin(2*np.pi*f*t)# Correct downsampling approachsample_rate=samples# Directly using samples as an integerdownsampling_factor=int(fs/sample_rate)# Downsampling factor# Downsample the signal by selecting every nth sampledownsampled_signal=signal[::downsampling_factor]# Recalculate L for the downsampled signal if neededL_downsampled=len(downsampled_signal)# Calculate the DFT of the downsampled signalsampled_signal_fft=fft(downsampled_signal)# Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signalsampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]n=len(downsampled_signal)T=1/sample_rate# Recalculate the sampling period with the selected sample rate# After calculating the magnitude of Fourier coefficientssampled_magnitude=np.abs(sampled_signal_fft)[:L_downsampled//2]# Create a new frequency vector that includes the replicated frequencies# First, calculate the original frequency bins for the positive frequenciesoriginal_sampled_freqs=fftfreq(L_downsampled,T)[:L_downsampled//2]shifted_freqs_negative=original_sampled_freqs-1/T# Shift for -1/Tsshifted_freqs_positive=original_sampled_freqs+1/T# Shift for +1/Ts# Concatenate the original and shifted (replicated) frequencies and magnitudes# This includes the original frequencies, and the replications at -1/Ts and +1/Tsreplicated_freqs=np.concatenate([shifted_freqs_negative,original_sampled_freqs,shifted_freqs_positive])replicated_magnitude=np.concatenate([sampled_magnitude,sampled_magnitude,sampled_magnitude])# Sort the replicated frequencies and magnitudes in ascending order for proper plottingsorted_indices=np.argsort(replicated_freqs)sorted_replicated_freqs=replicated_freqs[sorted_indices]sorted_replicated_magnitude=replicated_magnitude[sorted_indices]# Generate an array of sample indices for the downsampled signalsample_indices=np.arange(len(downsampled_signal))# Normalize the frequency valuesnormalized_freqs=original_sampled_freqs/(1/T)# Since you previously concatenated and sorted for replication, ensure to apply normalization there as wellnormalized_replicated_freqs=np.concatenate([(shifted_freqs_negative/(1/T)),# Normalize -1/Ts shifted frequenciesnormalized_freqs,# Already normalized original frequencies(shifted_freqs_positive/(1/T))# Normalize +1/Ts shifted frequencies])# Sort the normalized and replicated frequencies for proper plottingsorted_indices=np.argsort(normalized_replicated_freqs)sorted_normalized_replicated_freqs=normalized_replicated_freqs[sorted_indices]sorted_replicated_magnitude=replicated_magnitude[sorted_indices]# Adjust the sampling interval for selecting frequencies and magnitudessampled_indices=np.arange(0,len(sorted_normalized_replicated_freqs),selected_N)sampled_normalized_replicated_freqs=sorted_normalized_replicated_freqs[sampled_indices]sampled_replicated_magnitude=sorted_replicated_magnitude[sampled_indices]# Create the plotsfig,axs=plt.subplots(1,2,figsize=(15,5))# Time-domain signal plotaxs[0].scatter(sample_indices,downsampled_signal,color='#00CC96',s=12)axs[0].set_title('Sampled Signal')axs[0].set_xlabel('Sample Number')axs[0].set_ylabel('Amplitude')axs[0].grid(True)axs[1].stem(sampled_normalized_replicated_freqs,sampled_replicated_magnitude,linefmt='#1F77B4',markerfmt='o',basefmt=" ")axs[1].set_title('Sampled Normalized Fourier Transformation')axs[1].set_xlabel('Normalized Frequency (f/fs)')axs[1].set_ylabel('Magnitude')axs[1].grid(True)plt.tight_layout()plt.show()# Define interactive widgetsfrequency_slider=widgets.IntRangeSlider(value=[5,40],min=1,max=51,step=1,description='Frequency range (Hz):',layout=Layout(width='90%'),style={'description_width':'initial'})samples_slider=widgets.IntSlider(value=500,min=100,max=1000,step=25,description='Samples:',layout=Layout(width='90%'),style={'description_width':'initial'})N_slider=widgets.IntSlider(value=1,min=1,max=25,step=1,description='N:',layout=Layout(width='90%'),style={'description_width':'initial'})html_label=widgets.HTML(value=""" <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Properties of the Discrete Fourier Transform in Digital Signal Processing</h2> """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,frequency_slider,samples_slider,N_slider],layout=vbox_layout)out=widgets.interactive_output(plot_signals,{'f_range':frequency_slider,'samples':samples_slider,'selected_N':N_slider})# Display the widgets and outputclear_output(wait=True)# Clear the previous plotdisplay(ui,out)
To calculate the energy or power of the waveform \(x(t)\), depending on the case of the signal, it holds
\[
E_{X}=\int_{-\infty}^{\infty} x^{2}(t) d t=\int_{-\infty}^{\infty}|X(f)|^{2} d f
\]
\[
P_{X}=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{-T / 2}^{T / 2} x^{2}(t) d t=\int_{-\infty}^{\infty} S_{X}(f) d f
\]
where for power signals \(S_Χ(f)\) is the power spectral density (PSD) of \(x(t)\). For discrete-time signals resulting from sampling of \(x(t)\) with period \(T_s\), the corresponding relations for calculating the energy or power become
A simple way to estimate the power spectral density of the waveform \(x(t)\) is to take the DTFT of the signal’s samples and then square the magnitude of the result. This estimator is called a periodogram. The periodogram of a finite-length \(L\) signal \(x[n]\) is defined as
where \(X_d(\phi)\) the DTFT of the signal. As the length \(L\) tends to infinity, the periodogram \(P_{xx}(f)\) approaches the power spectral density \(S_Χ(f)\). The calculation of the periodogram at a finite number of frequencies \(kf_s/N\), \(k=0\), \(1\), \(…\), \(N\) gives
importnumpyasnpimportmatplotlib.pyplotasplt# ==============================================================================# 2.1 Create a Sinusoidal Signal# ==============================================================================Fs=2000# Sampling frequency in HzTs=1/Fs# Sampling period in secondsT=0.1# Signal duration in secondst=np.arange(0,T,Ts)# Time vector for signalA=1# Signal amplitudex=A*np.sin(2*np.pi*100*t)# Generate sinusoidal signalL=len(x)# Length of the signal# Plot the sinusoidal signal in time domainplt.figure()plt.plot(t,x)plt.title('Sinusoidal Signal')plt.xlabel('Time (s)')plt.ylabel('Amplitude')plt.show(block=False)plt.pause(1)# Pause to view the plot, press any key in the console to continue
% ==============================================================================% 2.1 Create a Sinusoidal Signal% ==============================================================================Fs=2000;% Sampling frequency in HzTs=1/Fs;% Sampling period in secondsT=0.1;% Signal duration in secondst=0:Ts:T-Ts;% Time vector for signalA=1;% Signal amplitudex=A*sin(2*pi*100*t);% Generate sinusoidal signalL=length(x);% Length of the signalplot(t,x);% Plot the sinusoidal signal in time domaintitle('Sinusoidal Signal');% Title for the plotxlabel('Time (s)');% X-axis labelylabel('Amplitude');% Y-axis labelpause;% Pause to view the plot, press any key to continue
Plot the discrete Fourier transform of the sine wave signal#
Code
# ==============================================================================# 2.2 Plot Fourier Transform (FT) of the Signal# ==============================================================================N=1*L# Length of Fourier TransformFo=Fs/N# Frequency resolutionFx=np.fft.fft(x,N)# Discrete Fourier Transform (DFT) of the signalfreq=np.arange(0,N)*Fo# Frequency vector# Plot the magnitude of the DFTplt.figure()plt.plot(freq,np.abs(Fx))plt.title('FFT of the Signal')plt.xlabel('Frequency (Hz)')plt.ylabel('Magnitude')plt.show(block=False)plt.pause(1)# Pause to view the plot, press any key in the console to continueplt.axis([0,100,0,L/2])plt.pause(1)
% ==============================================================================% 2.2 Plot Fourier Transform (FT) of the Signal% ==============================================================================N=1*L;% Length of Fourier TransformFo=Fs/N;% Frequency resolutionFx=fft(x,N);% Discrete Fourier Transform (DFT) of the signalfreq=(0:N-1)*Fo;% Frequency vectorplot(freq,abs(Fx));% Plot the magnitude of the DFTtitle('FFT of the Signal');% Title for the plotxlabel('Frequency (Hz)');% X-axis labelylabel('Magnitude');% Y-axis labelpause;% Pause to view the plot, press any key to continueaxis([01000L/2]);% Set axis limitspause;% Pause, press any key to continue
% ==============================================================================% 2.3 Plot Signal Periodogram% ==============================================================================power=Fx.*conj(Fx)/Fs/L;% Calculate spectral densityplot(freq,power);% Plot the periodogramtitle('{\bf Periodogram}');% Title for the plot with bold charactersxlabel('Frequency (Hz)');% X-axis labelylabel('Power');% Y-axis label
# ==============================================================================# 2.4 Calculate Signal Power# ==============================================================================power_theory=A**2/2# Theoretical power based on signal amplitudedB=10*np.log10(power_theory)# Convert power to decibels (dB)power_time_domain=np.sum(np.abs(x)**2)/L# Calculate power in time domainpower_frequency_domain=np.sum(power)*Fo# Calculate power in frequency domain# Display calculated power valuesprint(f'Power (Theory): {power_theory}')print(f'Power (dB): {dB}')print(f'Power (Time Domain): {power_time_domain}')print(f'Power (Frequency Domain): {power_frequency_domain}')
% ==============================================================================% 2.4 Calculate Signal Power% ==============================================================================power_theory=A^2/2;% Theoretical power based on signal amplitudedB=10*log10(power_theory);% Convert power to decibels (dB)power_time_domain=sum(abs(x).^2)/L;% Calculate power in time domainpower_frequency_domain=sum(power)*Fo;% Calculate power in frequency domain% Display calculated power valuesdisp(['Power (Theory): ',num2str(power_theory)]);disp(['Power (dB): ',num2str(dB)]);disp(['Power (Time Domain): ',num2str(power_time_domain)]);disp(['Power (Frequency Domain): ',num2str(power_frequency_domain)]);
Next, you will apply what you have learned in a more complex example of signal generation that includes modulation and noise addition. For convenience, incomplete PYTHON/MATLAB code with comments is provided, which you need to complete and save.
Code
# ==============================================================================# Part 1: Create the signal# ==============================================================================importnumpyasnpimportmatplotlib.pyplotasplt# Close all figure windows______# Clear workspace# In Python, variables need to be explicitly deleted if needed. Typically, just overwrite or ignore.# Clear command window# This is not applicable in Python in the same way as in MATLAB.Fs=4000# Sampling frequency in HzTs=___# Sampling periodL=3000# Signal length (number of samples)T=L*Ts# Signal durationt=np.arange(0,L)*Ts# Timestamps of signal calculation# Signal creationx=np.sin(2*np.pi*200*t) \
+0.3*np.sin(2*np.pi*300*(t-2)) \
+np.sin(2*np.pi*400*t)# Plot signal x in timeplt.figure(1)plt.plot(t,x)plt.title('Time domain plot of x')plt.xlabel('t (sec)')plt.ylabel('Amplitude')plt.show(block=False)# block=False allows the script to continue while the plot is shownplt.pause(1)# Pause for 1 second; adjust as neededplt.axis([0,0.3,-2,2])# Adjust axesplt.pause(1)# Pause for 1 second; adjust as needed# Calculate Fourier TransformN=2**np.ceil(np.log2(L)).astype(int)# FFT lengthFo=___# Frequency analysis stepf=np.arange(0,N)*Fo# Frequency vectorX=___# DFT for N points, fill in with the correct function# Plot signal in frequency domainplt.figure(2)plt.plot(f[:N//2],np.abs(X[:N//2]))# Plot signal only for positive frequency valuesplt.title('Frequency domain plot of x')plt.xlabel('f (Hz)')plt.ylabel('Amplitude')plt.show(block=False)plt.pause(1)# For a two-sided signal plotplt.figure(3)f=f-Fs/2# Shift frequency vectorX=np.fft.fftshift(X)# Shift the zero to the centerplt.plot(f,np.abs(X))plt.title('Two sided spectrum of x')plt.xlabel('f (Hz)')plt.ylabel('Amplitude')plt.show(block=False)plt.pause(1)# Calculate power of signalpower=(X*np.conj(X))/N/Lplt.figure(4)plt.plot(f,power)plt.xlabel('Frequency (Hz)')plt.ylabel('Power')plt.title('Periodogram')plt.show(block=False)plt.pause(1)
% ==============================================================================% Part 1 Create the signal% ==============================================================================_________%closeallfigurewindows_________%clearworkspace_________%clearcommandwindowFs=4000;% sampling frequency 4000 HzTs=____;% sampling periodL=3000;% signal length (number of samples)T=L*Ts;% signal durationt=0:Ts:(L-1)*Ts;% timestamps of signal calculationx=sin(2*pi*200*t)... % sinusoidal signal of frequency 200 Hz+0.3*sin(2*pi*300*(t-2))... % plus sinusoidal signal of frequency 300 Hz+sin(2*pi*400*t);% plus sinusoidal signal of frequency 400 Hz% Plot signal x in timefigure(1)% open a new figure windowplot(t,x)% plot the signaltitle('Time domain plot of x')% plot titlexlabel('t (sec)')% x axis labelylabel('Amplitude')% y axis labelpause% pause, press any key to continueaxis([00.3-22])% x axis values from 0 to 0.3 sec and% y axis values from -2 to 2pause% pause, press any key to continue% Calculate Fourier TransformN=2^nextpow2(L);% FFT length% nextpow2 calculates next power of 2% grater or equal to LFo=____;% frequency analysisf=(0:N-1)*Fo;% frequency vectorX=________;% DFT for Ν points% Plot signal in frequency domain% As signal is real you can plot the one sided signal(positive frequency values)figure(2)% open a new figure windowplot(f(1:_____)),abs(X(1:_____)))% plot signal only for positive%frequency valuestitle('Frequency domain plot of x')% plot titlexlabel('f (Hz)')% set x axis labelylabel('Amplitude')% set y axis labelpause% pause, press any key to continue% in order to plot the 2-sided signal you can use fftshift matlab% function, type help fftshift for more detailsfigure(3)% open a new figure windowf=f-Fs/2;% shift frequency vector to the left by –Fs/2X=fftshift(X);% create 2-sided signalplot(f,abs(X));title('Two sided spectrum of x');xlabel('f (Hz)');ylabel('Amplitude')pause% pause, press any key to continue% Calculate power of signalpower=X.*conj(X)/N/L;% Calculate power of signalfigure(4)% open a new figure windowplot(f,power)% plot signal power to frequencyxlabel('Frequency (Hz)')% x axis labelylabel('Power')% y axis labeltitle('{\bf Periodogram}')% plot titlepause
# Callback function to update graphsdefupdate_graph(selected_frequencies_part1,selected_Fo):# Part 1: Create the signalFs=selected_frequencies_part1# Sampling frequency 1000 HzTs=1/Fs# Sampling periodL=1000# Length of signal (number of samples)T=L*Ts# Duration of signalt=np.arange(0,(L-1)*Ts,Ts)# Time vectorglobalnew_xnew_x=np.sin(2*np.pi*(selected_Fo-30)*t) \
+0.8*np.sin(2*np.pi*(selected_Fo+20)*(t-2)) \
+np.sin(2*np.pi*selected_Fo*t)# 60 Hz component# Plottingfig,axs=plt.subplots(4,1,figsize=(12,20))# Time domain plotaxs[0].plot(t,new_x,color='#00CC96')axs[0].set_title('Time domain plot of x')axs[0].set_xlabel('t (sec)')axs[0].set_ylabel('Amplitude')axs[0].grid(True)# Fourier transformdefnextpow2(i):n=1whilen<i:n*=2returnnN=nextpow2(L)# Length of Fourier transformFo=Fs/N# Frequency resolutionf=np.arange(0,N)*Fo# Frequency vectorX=np.fft.fft(new_x,N)# Compute DFT for N points# Frequency domain plotaxs[1].plot(f[1:N],abs(X[1:N]),color='#1F77B4')axs[1].set_title('Frequency domain plot of x')axs[1].set_xlabel('f (Hz)')axs[1].set_ylabel('Amplitude')axs[1].grid(True)# Shift frequencies to centerf=f-Fs/2X=np.fft.fftshift(X)# Two-sided spectrum of xf_shifted=f# Two-sided spectrum plotaxs[2].plot(f_shifted,abs(X),color='#1F77B4')axs[2].set_title('Two sided spectrum of x')axs[2].set_xlabel('f (Hz)')axs[2].set_ylabel('Amplitude')axs[2].grid(True)# Calculate powerpower=np.multiply(X,np.conj(X))/N/L# Periodogram plotaxs[3].plot(f_shifted,power.real,color='#1F77B4')axs[3].set_title('Periodogram')axs[3].set_xlabel('Frequency (Hz)')axs[3].set_ylabel('Power')axs[3].grid(True)plt.tight_layout()plt.show()# Create interactive widgetssingle_frequency_slider=widgets.IntSlider(min=100,max=2000,step=100,value=1000,description='Sampling Frequency (Fs):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)Fo_slider=widgets.IntSlider(min=40,max=400,step=10,value=100,description='Frequency (Fo):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,single_frequency_slider,Fo_slider],layout=vbox_layout)out=widgets.interactive_output(update_graph,{'selected_frequencies_part1':single_frequency_slider,'selected_Fo':Fo_slider})# Display the widgets and outputclear_output(wait=True)# Clear the previous plotdisplay(ui,out)
Complete the code to create the noise signal n using the randn function. The noise vector n should be the same size as the sine wave x from the first part. Plot the noise signal in the interval from 0 to 0.2 sec and scale from -2 to 2. Calculate the periodogram of n and plot the power spectral density of the noise signal. Add the noise signal and x to get the noisy signal s. Plot the noisy signal s in the time domain in the area from 0 to 0.2 sec and scale from -2 to 2 as well as its bilateral spectrum.
Fs=1000# συχνότητα δειγματοληψίας 1000 HzTs=1/Fs# περίοδος δειγματοληψίαςL=1000# μήκος σήματος (αριθμός δειγμάτων)T=L*Ts# διάρκεια σήματοςt=np.arange(0,(L-1)*Ts,Ts)# χρονικές στιγμές υπολογισμού του σήματοςdefnextpow2(i):# Compute the next highest power of 2n=1whilen<i:n*=2returnn# Function to update plotsdefupdate_plots(selected_frequencies_part2):new_x=np.sin(2*np.pi*(selected_frequencies_part2-30)*t)+0.8*np.sin(2*np.pi*(selected_frequencies_part2+20)*(t-2))+np.sin(2*np.pi*(selected_frequencies_part2)*t);rand_n=np.random.randn(np.size(new_x))# Plottingfig,axs=plt.subplots(4,1,figsize=(12,20))# Time domain plot of naxs[0].plot(t,rand_n,color='#00CC96')axs[0].set_title('Time domain plot of n')axs[0].set_xlabel('t (sec)')axs[0].set_ylabel('Amplitude')axs[0].grid(True)# Correction for N calculation using bitwise operatorN=2^nextpow2(L)Fo=Fs/Nf=(np.arange(0,N))*Fof_shifted=f-Fs/2rand_N=np.fft.fft(rand_n,N)rand_N=np.fft.fftshift(rand_N)power_n=np.multiply(rand_N,np.conj(rand_N))/N/L# Frequency domain plot of xaxs[1].plot(f_shifted,power_n.real,color='#1F77B4')axs[1].set_title('Frequency domain plot of x')axs[1].set_xlabel('f (Hz)')axs[1].set_ylabel('Amplitude')axs[1].grid(True)# Two sided spectrum of xs=new_x+rand_naxs[2].plot(t,s,color='#00CC96')axs[2].set_title('Two sided spectrum of x')axs[2].set_xlabel('t (sec)')axs[2].set_ylabel('Amplitude')axs[2].grid(True)# Two sided spectrum of sS=np.fft.fft(s,N)S=np.fft.fftshift(S)axs[3].plot(f_shifted,np.abs(S),color='#1F77B4')axs[3].set_title('Two sided spectrum of s')axs[3].set_xlabel('f (Hz)')axs[3].set_ylabel('Magnitude')axs[3].grid(True)plt.tight_layout()plt.show()# Create the slider widgetfs_slider=widgets.IntSlider(value=1000,min=100,max=2000,step=100,description='Sampling Frequency (Hz):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,fs_slider],layout=vbox_layout)out=widgets.interactive_output(update_plots,{'selected_frequencies_part2':fs_slider})# Display the sliderclear_output(wait=True)# Clear the previous plotdisplay(ui,out)
Complete the code for creating a sine wave signal of 100 Hz frequency and multiply it with the previous signal s. The two signals should be of the same size. Plot the result in the time domain in the area from 0 to 0.2 sec and scale from -2 to 2 as well as in the frequency domain using the fftshift function.
# Part 3. Πολλαπλασιασμός σημάτων# Συμπληρώστε τον κώδικα δημιουργίας ενός ημιτονοειδούς σήματος συχνότητας# 100 Hz και πολλαπλασιάστε με το προηγούμενο σήμα s.# Τα δύο σήματα θα πρέπει να είναι του ίδιου μεγέθους.# Σχεδιάστε το αποτέλεσμα στο πεδίο του χρόνου στην περιοχή 0 έως 0.2 sec# και κλίμακα από -2 έως 2 καθώς και στο πεδίο της συχνότητας# χρησιμοποιώντας τη συνάρτηση fftshift.defupdate_plots(selected_frequencies_part3):Fο=selected_frequencies_part3z=np.sin(2*np.pi*Fο*t)new_x=np.sin(2*np.pi*30*t) \
+0.8*np.sin(2*np.pi*80*(t-2)) \
+np.sin(2*np.pi*60*t)# 60 Hz componentL=1000# Length of signalrand_n=np.random.normal(0,1,np.size(new_x))# Example random noises=new_x+rand_ny=np.multiply(z,s)# Plottingfig,axs=plt.subplots(2,1,figsize=(12,20))# Time domain plotaxs[0].plot(t,y,color='#00CC96')axs[0].set_title('Time domain plot of y')axs[0].set_xlabel('t (sec)')axs[0].set_ylabel('Amplitude')axs[0].grid(True)# Fourier transformdefnextpow2(i):n=1whilen<i:n*=2returnnN=2^nextpow2(L)Fo=Fs/Nf=(np.arange(0,N))*FoY=np.fft.fft(y,N)f=f-Fs/2Y=np.fft.fftshift(Y)axs[1].plot(f,np.abs(Y),color='#1F77B4')axs[1].set_title('Frequency domain plot of y')axs[1].set_xlabel('f (Hz)')axs[1].set_ylabel('Amplitude')axs[1].grid(True)plt.tight_layout()plt.show()# Create the slider widget for FoFo_slider_step3=widgets.IntSlider(value=100,min=10,max=200,step=10,description='Fo (Hz):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,Fo_slider_step3],layout=vbox_layout)out=widgets.interactive_output(update_plots,{'selected_frequencies_part3':Fo_slider_step3})# Display the sliderclear_output(wait=True)# Clear the previous plotdisplay(ui,out)
Write a Python spectral analysis function, similar to signal.welch(): it will accept as input a real signal vector as well as the sampling frequency, \(F_s\), and will plot the one-sided spectral density of the signal in the range \([0-F_s/2)\). The signal will be segmented into sections of length equal to the power of \(2\) closest to \(1/8\) of its total length, but not less than 256. The sections will overlap by \(50%\). The last section, if shorter than the others, will be ignored. The spectrum of each section will be calculated with FFT and the mean value of all sections will be taken. The function should be tested with the signal from example 1.1 and the result compared with that of signal.welch().
Code
# ==============================================================================# Custom nextpow2 and pwelch functions# ==============================================================================defnextpow2(i):n=1whilen<i:n*=2returnndefpwelch(x,Fs):Ts=1/FsL=np.size(x)+1T=L*TsN=2^nextpow2(L)Fo=Fs/Nf=np.arange(0,N)*Fowindow_size=nextpow2(np.size(x)/8)if(window_size<256):window_size=256windows=np.size(x)//(window_size//2)-1indexer=np.arange(window_size)[None,:]+(window_size//2)*np.arange(windows)[:,None]windowed_x=x[indexer]avg_pwr=0forwindowinwindowed_x:window=window*np.hanning(np.size(window))L=np.size(window)+1T=L*TsN=2^nextpow2(L)Fo=Fs/Nf=np.arange(0,N)*Fowindow_fft=np.fft.fft(window,N)power=np.multiply(window_fft,np.conj(window_fft))/N/Lavg_pwr=avg_pwr+poweravg_pwr=avg_pwr/windowsfig,ax=plt.subplots()fig.set_size_inches(18.5,10.5)ax.plot(f[np.arange(0,N//2)],avg_pwr[np.arange(0,N//2)])ax.set(xlabel='Frequency (Hz)',ylabel='Power',title='Periodogram pwelch()')ax.grid()plt.show()returnf[np.arange(0,N//2)],avg_pwr[np.arange(0,N//2)]# ==============================================================================# Test the functions & plot the results# ============================================================================== Fs=______f1,Pxx1=pwelch(x,Fs)f2,Pxx2=signal.welch(x,fs=Fs)fig,ax=plt.subplots()fig.set_size_inches(18.5,10.5)ax.plot(f2,Pxx2)ax.set(xlabel='Frequency (Hz)',ylabel='Power',title='Periodogram signal.welch()')ax.grid()plt.show()
% ==============================================================================% Custom nextpow2 and pwelch functions% ==============================================================================functionN=nextpow2(i)n=1;whilen<in=n*2;endN=log2(n);% This returns the exponent N such that 2^N = n, aligning with MATLAB's built-in nextpow2 behavior.endfunction[f, avgPwr] = pwelch(x, Fs)Ts=1/Fs;L=length(x)+1;T=L*Ts;N=2^nextpow2(L);Fo=Fs/N;f=(0:N-1)*Fo;windowSize=nextpow2(length(x)/8);if(windowSize<256)windowSize=256;endwindows=floor(length(x)/(windowSize/2))-1;indexer=bsxfun(@plus,(1:windowSize)',(windowSize/2)*(0:windows-1));windowedX=x(indexer);avgPwr=0;fori=1:size(windowedX,2)window=windowedX(:,i).*hanning(size(windowedX,1));L=length(window)+1;T=L*Ts;N=2^nextpow2(L);Fo=Fs/N;f=(0:N-1)*Fo;windowFFT=fft(window,N);power=(windowFFT.*conj(windowFFT))/(N*L);avgPwr=avgPwr+power;endavgPwr=avgPwr/windows;figure;set(gcf,'Position',[100,100,1850,1050]);plot(f(1:floor(N/2)),avgPwr(1:floor(N/2)));xlabel('Frequency (Hz)');ylabel('Power');title('Periodogram pwelch()');gridon;f=f(1:floor(N/2));avgPwr=avgPwr(1:floor(N/2));end% ==============================================================================% Test the functions & plot the results% ==============================================================================Fs=_______;% Define your sampling frequencyx=randn(1,10000);% Example data, replace with your actual data[f1,Pxx1]=pwelch(x,Fs);[f2,Pxx2]=pwelch(x,Fs,[],[],Fs);% MATLAB's built-in functionfigure;set(gcf,'Position',[100,100,1850,1050]);plot(f2,Pxx2);xlabel('Frequency (Hz)');ylabel('Power');title('Periodogram signal.welch()');gridon;
defnextpow2(i):n=1whilen<i:n*=2returnndefpwelch(x,Fs):Ts=1/FsL=np.size(x)+1T=L*TsN=2^nextpow2(L)Fo=Fs/Nf=np.arange(0,N)*Fowindow_size=nextpow2(np.size(x)/8)if(window_size<256):window_size=256windows=np.size(x)//(window_size//2)-1indexer=np.arange(window_size)[None,:]+(window_size//2)*np.arange(windows)[:,None]windowed_x=x[indexer]avg_pwr=0forwindowinwindowed_x:window=window*np.hanning(np.size(window))L=np.size(window)+1T=L*TsN=2^nextpow2(L)Fo=Fs/Nf=np.arange(0,N)*Fowindow_fft=np.fft.fft(window,N)power=np.multiply(window_fft,np.conj(window_fft))/N/Lavg_pwr=avg_pwr+poweravg_pwr=avg_pwr/windowsreturnf[np.arange(0,N//2)],avg_pwr[np.arange(0,N//2)]# Function to update plots based on slider valuedefupdate_plots(Fs):T=1/Fs# Update sampling periodt1=np.arange(0,L)*T# Update time vector# Recompute signal x with new sampling frequencylast_x=np.sin(2*np.pi*30*t1)+0.8*np.sin(2*np.pi*80*(t1-2))+np.sin(2*np.pi*60*t1)# Compute pwelchf1,Pxx1=pwelch(last_x,Fs)# Compute signal.welchf2,Pxx2=signal.welch(last_x,fs=Fs)# Plotfig,axs=plt.subplots(2,1,figsize=(18.5,20))# Plot custom pwelchaxs[0].plot(f1,Pxx1)axs[0].set(xlabel='Frequency (Hz)',ylabel='Power',title='Periodogram pwelch()')axs[0].grid()# Plot signal.welchaxs[1].plot(f2,Pxx2)axs[1].set(xlabel='Frequency (Hz)',ylabel='Power',title='Periodogram signal.welch()')axs[1].grid()plt.tight_layout()# Create slider for FsFs_slider=widgets.IntSlider(value=500,min=100,max=2000,step=100,description='Sampling Frequency (Fs):',layout=Layout(width='90%'),style={'description_width':'initial'},continuous_update=False)html_label=widgets.HTML(value=""" """)# Display the sliders and outputvbox_layout=Layout(display='flex',flex_flow='column',align_items='center')ui=widgets.VBox([html_label,Fs_slider],layout=vbox_layout)out=widgets.interactive_output(update_plots,{'Fs':Fs_slider})# Display the sliderclear_output(wait=True)# Clear the previous plotdisplay(ui,out)